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Abstract 

Free energy of crystal phases is commonly evaluated by thermodynamic integration (TDI) along 
a reversible path that involves an external potential. A persistent problem in this method is 
that a significant hysteresis is observed due to differences in the center of mass position of the 
crystal phase in the presence and absence of the external potential. To alleviate this hysteresis, a 
constraint on the translational degrees of freedom of the crystal phase is imposed along the path 
and subsequently a correction term is added to the free energy to account for such a constraint. 
In this work, we propose a new methodology termed as error-biased Bennett Acceptance ratio 
(EBAR) method that effectively solves this problem without the need to impose any constraint. 
This method is simple to implement and it does not require any modification to the path. We show 
the applicability of this method in the computation of crystal-melt interfacial energy by cleaving 
wall method [J. Chem. Phys., 118, 7651 (2003)] and bulk crystal-melt free energy difference by 
constrained fiuid A-integration method [J. Chem. Phys., 120, 2122 (2004)] for a model potential 
of silicon. 
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I. INTRODUCTION 



Accurate evaluation of free energy of a crystal phase by simulation is important to pre- 
dict many important thermodynamic properties including the melting temperature,- the 
relative stability of different crystal phases^ (i.e., polymorphism), and the interfacial free 
energy of crystal phases. ^'^ A common method of free energy evaluation is thermodynamic 
integration (TDI), in which an external potential is imposed on the crystal phase to en- 
sure thermodynamic reversibility,- and the method requires that the resulting free energy 
change be computed accurately. However, the external potential causes center of mass (CM) 
position of the crystal phase to change relative to that in the absence of the external po- 
tential and this leads to hysteresis in the free energy computation. This problem has been 
reported to occur in many TDI methods including Einstein crystal method,-"^ constrained 
fluid A-integration method,-^'! surface free energy calculation of crystal phases,- and direct 
computation of crystal-melt interfacial energy by cleaving wall method.- 

A variety of constraints are imposed on the translational degrees of freedom of the crys- 
tal phase to address this problem. In Einstein crystal method, CM is fixed during ther- 
modynamic integration along the path.- In direct computation of crystal-melt free energy 
difference, CM is either fixed^ or confined to a small region.- A fixed particle constraint was 
also proposed for this method.'' In surface free energy calculation of Au[110] crystal phase,- 
CM velocity was artificially controlled by the (otherwise) dormant external potential during 
Molecular Dynamics (MD) simulation. In the computation of crystal-melt interfacial en- 
ergy 7, two crystal layers were immobilized by assigning them an infinite mass throughout 
the four stage integration path.- It can be seen from the above examples that although 
cause of the hysteresis is the same, every TDI path requires its own ad-hoc procedure to 
constrain the translational degrees of freedom of the crystal phase. Further, the use of a 
constraint entails computation of the correction term to obtain free energy difference be- 
tween the unconstrained phases. In case of Einstein crystal method,- the correction term 
can be computed analytically and the effort required is minimal. When CM is confined 
to a small region in constrained fluid A-integration method, the correction term needs to 
be computed by separate simulations.- In cleaving wall method,- 7 was computed for two 
different system sizes in order to determine the effect of the constraint. Thus, the estimation 
of the correction term is, in itself, computationally expensive in many cases. 
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In this article, we propose a new methodology termed as error-biased Bennett acceptance 
ratio (EBAR) method (to be explained in sec. II) to efficiently compute the free energy 
change resulting from imposition of the external potential on the crystal phase, without re- 
quiring the use of any constraint. We demonstrate two applications of this methodology in 
sec. III. We ffist compute the free energy of the crystal phase due to repulsive wall potential 
in cleaving wall method.-!^ Next, we compute the free energy difference resulting from the 
imposition of an attractive Gaussian well potential on the crystal phase in constrained fluid 
A-integration method.-*^ Both calculations are performed for the Stillinger and Weber- po- 
tential of silicon. Finally, we summarize our results and discuss further possible applications 
of this methodology in sec. IV. 



II. ERROR-BIASED BENNETT ACCEPTANCE RATIO METHOD 

In this section, we ffist describe the Bennett acceptance ratio methoci^ (BAR) which 
is routinely employed to compute the free energy difference between two states.— >i2i>i^ We 
then describe the error biased Bennett acceptance ratio (EBAR) method in the context 
of the present problem. According to the BAR method,— Helmholtz free energy difference 
AF = Fi — Fq between two equilibrium states '0' and '1' is given by the following equation: — 

AF = log — — — + C-log— , (1) 

Eo/(0i - 0o) - <-) no 

where C is a constant, f{x) = 1/(1 + e^) is the fermi function, J2o J2i represent the 
sums over fermi functions sampled in '0' and '1' ensembles, respectively. The symbols 0o 
and 01 in Eq. ([1]) represent the total configurational energies and hq and ni are the number 
of sampled fermi functions in the above two ensembles. Please note that we have expressed 
AF, 00, 01, and C in units of fc^T and we will follow this convention throughout unless 
otherwise stated explicitly. The error (in units of ksT) in the free energy estimate is given 
by 0"^ = ((AF — AA)"^), where AA is the expectation (correct) value of the free energy 
difference. When we are in the large sample regime (i.e., both J2o and J2i are reasonably 
accurate), the variance can be approximated 

2 {p)o - {f)i , - m ... 

where (/)o is the ensemble average of /(0i — 0o — C) evaluated in ensemble and (/)i is 
the ensemble average of /(0o — 0i + C) evaluated in ensemble 1. Bennett showed that the 



value of C = Cb which minimizes is given by the following expressions^ 

AF = CB-log^. (3) 

no 

This value of C corresponds to the condition^° that dcr'^/dC = 0. Thus, in the BAR method, 
the optimum estimate of AF is obtained by solving Eqs. ([1]) and ([3]) simultaneously. The 
condition given by Eq. ([3]) can also be expressed as J2o = The important requirement 
for the applicability of the BAR method is that the large sample regime should be achieved 
in both the ensembles. 

Considering the present problem, let's denote the crystal phase with and without an 
external potential as state 1 and 0, respectively. In such a case, the instantaneous config- 
urational energies in the two states are given by (po=U and 4>i=U + Uext, where U is the 
potential energy due to interactions between the particles and Uext is the external potential 
energy. Because of the influence of Uext, average position of center of mass of the crystal 
phase will be different in the two states. In order to sample the perturbation energy 0o ~ 4>i 
efficiently in state 1 the following two conditions are necessary, (i) The important configu- 
rations must overlap to a large extent with those of state 0, which is possible if the external 
potential does not affect the crystal structure significantly, (ii) Further, the external poten- 
tial must be sufficiently strong so that the CM is localized at the average position Ri. If 
this latter condition is not satisfied, the number of important configurations (corresponding 
to all possible CM positions) becomes too large and cannot be sampled in a finite length 
simulation. In state 0, the important values of the perturbation energy 0i — 0o are those 
for which the CM position is Ri (corresponding to the state 1) and are therefore not likely 
to be sampled in a simulation of reasonable length because the CM position can fiuctuate 
wildly in state due to the absence of the external potential. This leads to a large error in 
the estimate of J2o Eq. ([T]) for a given value of C and uq. 

When conditions (i) and (ii) are satisfied, the value of J2i in Eq. ([1]) can be estimated 
accurately in a simulation of reasonable length. However, due to poor estimate of J2o, the 
large sample regime is not achieved aX C = Cb and the BAR method fails. To circumvent 
this difficulty, we propose the EBAR method, in which we choose an appropriate value of 
C = Ce (away from Cb) such that the estimate of both J2o and J2i are reasonably accurate 
and hence the large sample regime is achieved. For this purpose, it is important to note the 
following two points: 
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When the value of C is chosen such that X^o ^ log So ~ log(So))~ where (Z^o) 



IG and 



is the expectation (accurate) value of X^o ^ given C and no (see Fig. 4 of Ref. 
the accompanying discussion). In such a case, the main source of error will be due to 
which we expect to be sufficiently accurate when the conditions (i) and (ii) given above are 
satisfied. 

(b) In order to locate the exact value of C where a large sample regime will occur, we first 
consider two limiting cases. When C approaches C+oo, such that x = (0o ~ + C) +00, 
f{x) ~ and /(— x) ~ 1. In this case, Eq. ([1]) reduces to the single stage free energy 
perturbation (FEP) formula AF = log(exp[— (00 — 0i)])i, where (■ ■ ■)! is the NVT ensemble 
average in state 1. In this limit da'^/dC = 0, however, we shall be necessarily in the small 
sample regime since we are not using any data from the state in the estimation of AF [the 
FEP formula in state 1 can be considered as a limiting caseAfi of the acceptance ratio formula 
in Eq. ([1]) as no 0]. Also at C = Cb, dcr^/dC = according to the BAR method, but the 
result corresponds to the small sample regime due to the poor estimate of Xo? as discussed 
before. For an intermediate value of C = Ce {Cb < Ce < C+oo), the magnitude of do^ jdC 
must attend a local maximum. The large sample regime will occur in the immediate vicinity 
of Ce-, since as we move away from Ce in either direction (C Cb or C ^ C+oo), the 
magnitude of da"^ /dC decreases to and we approach the small sample regime. 

To summarize the EBAR method, we choose C = Ce such that the magnitude of 
\da'^/dC\ is at its maximum and Xo ^ 1- It may be noted that a local maximum of 
\da'^ /dC\ will also exist when Xi ^ 1, but it will not correspond to the large sample regime 
due to poor estimate of Xo- In order that the EBAR method be successful, the estimation 
of Xi niust be sufficiently accurate. This is ensured by the conditions (i) and (ii) mentioned 



above. As in Ref. 
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the above methodology can be extended straightforwardly to isother- 
mal isobaric (NPT) ensemble. In this case, the values of the perturbation energies in above 
equations will be sampled during NPT simulations and Eqs. ([I])-([3]) yield the Gibbs free 
energy difference AG between states and 1 instead of the Hemholtz free energy difference 
AF. 
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III. APPLICATIONS 



In what follows, we demonstrate two applications of the EBAR method for the Stillinger- 
Weber (SW) potential- of silicon at the previously reported melting point, i.e., at T* = 
0.0667 (1678 K) and P* = 0^ (quantities with the superscript * are dimensionless) . 

A. Cleaving wall method 

In cleaving wall method, bulk crystal and melt phases are combined reversibly to form an 
interface and the work required for this change per unit interfacial area yields the crystal-melt 
interfacial energy 7.- This technique is sufficiently precise so as to resolve the anisotropy of 7 
and has been applied to study anisotropy of interfacial energies of hard-spheres,— soft-sphere 
potential, and to model potentials of silicon^ and water.— In the first stage of this process, 
crystal phase is cleaved at a predetermined plane so that no particle can cross that plane. 
The cleaving is done by introducing a repulsive wall consisting of one or more ideal crystal 
layers. These layers are chosen to have the same orientation as that of the interface being 
studied. This process is prone to hysteresis due to center of mass motion of the crystal 
phase since it changes the relative distance between the wall and the crystal layers. To 
prevent the hysteresis, Davidchack and Laird^ immobilized particles in the two layers of the 
crystal phase sufficiently far from the cleaving plane while applying the cleaving potential. 
To study the effect of this fixed layer constraint, thermodynamic integration was performed 
with the same interfacial area but with fewer crystal layers and the resulting value of 7 was 
found to be in agreement within error bars with that of the larger system, which indicated 
that no correction term was necessary.- However, computation of 7 for two different system 
sizes is computationally expensive since it involves implementation of all four stages of the 
cleaving wall method. We, therefore, explored the possibility that using EBAR method, we 
may be able to compute free energy accurately without the need to immobilize particles in 
the crystal layers. 

We performed the computation for (111) orientation of the silicon crystal phase, since 
for this orientation the layers on the opposite sides of the cleaving plane are not symmetric. 
As a result, when the CM of the crystal phase fluctuates, it results in fluctuation of the 
cleaving potential, thus leading to hysteresis. For (100) orientation, this problem does not 
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occur because of the symmetry of layers on the opposite sides of the cleaving plane. As in 
Ref. y, we employed crystal phase with 3024 particles and a simulation box with dimensions 
of = 3-\/2a, Ly = 3v^L5a, and = lAy/Sa, where a = {S/pcY^'^ is the unit cell length 
of the crystal phase and pc = 0.452(T~^ is the crystal density. The crystal phase contained 
84 layers of (111) orientation in the z-direction, with 36 particles in each layer. The total 
configurational energy at a given value of z is given by 

0(r^z) = f/(r^)+[/Ur^z), (4) 

where = (ri,---,r7v) is the instantaneous configuration of the particles, Uext is the 
cleaving potential exerted by the wall, and z is the distance between the wall and the cleaving 
plane. The cleaving potential consisted of a repulsive two-body term and a 3-body term 
derived from Stillinger and Weber potential as explained in Ref. [8. The cleaving wall was 
constructed of two ideal crystal layers of (111) orientation with 36 particles in each layer .- 
The cleaving planes were located at two boundaries of the simulation box in z-direction. 
The Helmholtz free energy change for cleaving of the crystal phase is given by the following 



expression :- 



^/ dF 

AF = I dz— (5) 

'■f , /dd\ 
dz 



, dz I 

where (■ ■ ■) is canonical ensemble average at a particular value of z. Following Ref. 
we chose the initial and final positions of the cleaving walls as Zi = 1.80 and Zf = 0.75, 
respectively. We performed canonical ensemble Monte Carlo (MC) simulations at various 
values of z with 50000 MC steps for equilibration and 2 x 10^ steps for production run. Each 
MC step consisted of 3024 trial displacement moves. The maximum value of the attempted 
displacement during the trial moves was adjusted during the equilibration period to have 
an acceptance ratio of nearly 50 %. The integrand in Eq. ([5]) and the perturbation energy 
[(01 — 0o) or (00 — 0i) required in Eq. ([1]) ] was sampled after every MC step. 

Figure [U shows the negative of the integrand in Eq. (5) per unit area (obtained without 
applying any constraint) as a function of the distance z. The results and the corresponding 
error bars reported in Fig. [1] and elsewhere represent the mean and the standard error of 
the estimates obtained from block averages during the production run. The ordinate in the 
plot is the magnitude of the force between the wall and the particles, which increases with 



decreasing z because of the repulsive nature of the walL The integrand shows large hysteresis 
throughout the path in the absence of the any constraint. However, when we apply the fixed 
layer constraint the hysteresis disappears as seen in the inset of Fig. [TJ When performing 
simulations with the constraint,^ particles in two central layers (layers numbering 42 and 43) 
in the z-direction were fixed, i.e., trial moves attempting to displace the particles in these 
layers were rejected with certainty during MC simulation. A qualitative comparison of the 
plots shown in Fig. [1] confirms that fluctuations in the CM leads to the hysteresis. Note that 
due to the effect of the constraint on the free energy, the values of the integrand in the inset 
of Fig. [Hare significantly different from those in the main plot. Due to this hysteresis, the 
BAR method is expected to fail and hence we have applied the EBAR method as explained 
below. 

In order to compute the free energy difference, we consider state as the crystal phase 
with a cleaving wall distances of 2;=!. 8 which is equal to the cut-off distance of the external 
potential. Thus, in state no external potential acts on the system. As we decrease the 
value of z, the influence of the external potential acting on the system increases. In Figs. [2] 
and [31 we report the relevant details of the EBAR and the BAR methods when using data 
from two simulations at z = 0.75 (state 1) and 1.8 (state 0). As seen in Fig. [2|, the BAR 
result, which corresponds to the condition da"^ /dC = (or equivalently J2o = is in the 
small sample regime,— since J2o = J2i « 1- Also, from the slope of the curve at C = Cb in 
Fig. [31 we find that dAF/dC = —1 which again confirms that the BAR result is in the small 
sample regime.— Thus the optimization scheme fails and the BAR method result is far off 
from the actual value, as seen in Fig. [31 In the EBAR method, we choose a value of Ce (see 
Fig. [21) such that the magnitude of da'^ /dC is maximum and J2o ^ 1- The corresponding 
value of AF (see Fig. [31) compares well with the accurate result reported in Ref. Is using the 
BAR method. At a large positive value of C > 220, the value of AF/A computed using 
the acceptance ratio formula Eq. approaches 0.0175 J/m^, which corresponds to the 
FEP formula in state 1. Thus, the EBAR method (see Fig. [31), yields more accurate results 
compared to single stage FEP method. This is expected since the EBAR method utilizes 
data from both the ensembles unlike the FEP method which relies on data from ensemble 1 
alone. 

Figure [H compares the computational effort required to obtain AF/A using different 
methodologies (see also Table [11). When considering more than two simulations {Ng > 2), 
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the EBAR method is apphcable only to the last interval and hence is combined with the 
BAR method to obtain total free energy difference. As an example, when considering 4 
simulations {N^ = 4) at z = 0.75, 0.8, 0.85, and 1.8, the result reported under the combined 
method (C) (see Table 1 and Fig. Hj) is a sum of BAR results from z = 0.75 to 0.85 and the 
EBAR result from z = 0.85 to 1.8. Figure H] shows that the EBAR method yields accurate 
value with just two simulations while the BAR method requires at least 9 simulations. It can 
also be seen that numerical integration using Gaussian Quadrature (GQ) technique using 3 
simulations is close to the correct value, however the size of the error bar is still relatively 
large compared to the EBAR result. This clearly indicates that GQ will require 4 or more 
simulations to achieve the desired accuracy. 

We also tested for the convergence of the EBAR method in Fig. H] by using simulation 
data at intermediate z values. For example for A^^ = 2, 4, and 6, the EBAR method was 
applied to the intervals from z = 0.75 to 1.8, z = 0.85 to 1.8, and z = 1.0 to 1.8, respectively. 
As the width of the interval reduces, the extent of the overlap between the configurations 
in state 1 and state will increase. However, the CM becomes less localized in state 1 due 
to reduced strength of the external potential at higher value of z and hence the estimation 
of J2i becomes less accurate for a given rii as discussed in Sec. II. Note also that the size of 
the error bars in Fig. [1] is larger for z > 0.75 compared to that for z = 0.75, indicating a 
larger hysteresis due to CM fluctuations. Since the accuracy of the EBAR method depends 
on X]i5 "we find that it becomes less accurate (although marginally) for smaller intervals (see 
EBAR/Combined results in Table [land Fig. 11). 

B. Constrained fluid A-integration 

Grochola^ introduced constrained fluid A-integration method to compute the free energy 
difference between the crystal and the melt phases. A major advantage of this method 
compared to traditional TDI methods^ is that it directly calculates the free energy difference 
between the melt and the crystal phases by constructing a reversible path between the two 
phases. Thus, there is no need to connect the crystal and the melt phases separately to 
reference phases of known free energy and this gives more flexibility in terms of designing 
the reversible path. Grochola's method has been applied to calculate melting temperature 
of complex potentials for Sodium Chloride,— Benzene and trizole.— Further, the method 
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was extended to isothermal isobaric ensemble,- and to computation of melting temperature 
of binary mixtures^fl and interfacial free energy of crystal phases.- 

Here we compute Gibbs free energy difference between crystal and melt phases of silicon 
by NPT version'' of the constrained fluid A— integration method. The expressions for the 
potential energy for the 3 stages of the reversible path are as given below: 

0i(Ai) = (l-r7Ai)?7, (6) 

02(A2) = (l-r/)t/ + A2t/e.t, (7) 

and 

03(A3) = [(1 -r]) + ^^r,]U + (1 - \-i)Ue.u (8) 

where U is the potential energy due to interactions between the system particles, ?7 is a 
parameter controlling the extent to which strength of interaction is reduced in the first 
stage.-*^ Ai, A2, and A3 are the parameters characterizing the three stages. The Gaussian 
external potential imposed during the second and the third stage is given by^'^ Uext = 
Z]j Z]fc '2 6xp(— for^^^), where the summation with respect to i is taken over all the particles, 
is the distance between the ith particle and fcth well, and the summation with respect 
to k is taken over all Gaussian potential wells within a certain cutoff distance of the zth 
particle. A constraint on the maximum possible volume Vm along the path is imposed^ to 
ensure thermodynamic reversibility. The Gibbs free energy change for the third stage of the 
path is given by:- 

AG3= f^dX3{vU-Ue.t), (9) 
Jo 

where (■ ■ ■) represents the isothermal-isobar ic ensemble average at a given value of A3. 
Similar expressions apply for AGi and AG2.- 

The simulation at a given state point along the path was performed in NPT ensemble 
with N = 1000 particles confined to a cubic simulation box under periodic boundary con- 
ditions. At a given state point, we used 15000 MC steps for equilibration and 20000 steps 
for production, except at the end of 3rd stage, where 10^ MC steps were used for produc- 
tion run from A3 = 0.94 to 1. In each MC step, we attempted, on average, two volume 
change moves and 1000 particle displacement moves. The size of the attempted changes 
was adjusted during equilibration steps so as to achieve an acceptance ratio of about 50%. 
During production run, we sampled the integrand for TDl method and the perturbation 
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energy for BAR method after every MC step. The Gaussian well parameters were chosen to 



a 



20 



The 



be a = — 1.892e and b = 8.0a in accordance with the criteria mentioned in Refs. 
parameter t] in Eqs. (jH])^® was assigned a value of 0.9 following earlier work.^''' Further, 
we chose the maximum volume to be ^ = A^/0.4 at T* = 0.0667 and P* = so that Vm 
does not affect the free energy of the crystal or the melt phases.- 

We found no hysteresis for the first two stages of the path as in earlier studies.^ In the third 
stage, however, we found hysteresis between A3 = 0.99 to 1.0 as seen in the inset of Fig. [51 
The reason for this hysteresis is that as (A3 — 1), the influence of the external potential on 
the crystal phase becomes negligible as can be seen in expression of 03 in Eq. (JHl), which 
results in fluctuations of the CM position. Note that in computations performed by Grochola 
(see Fig. 6 of Ref. O), no hysteresis was seen in the third stage, because zero CM velocity 
was maintained during MD simulations. In order to compute the free energy difference in 
stage 3, we denote state corresponding to A3 = 1 as state 1, in which no external potential 
acts on the crystal phase. Note that the role of states and 1 is reversed compared what 
we discussed earlier since the external potential is acting in state (A3 < 1) in the present 
case. As A3 is decreased, the influence of the external potential increases. We compared 
BAR, EBAR, GQ methods on the basis of number of simulations (A^^) performed at various 
values of A3 between 0.9 and 1 to obtain a given result. For Ns > 2, the EBAR result is 
applicable only to the last interval and hence we combine it with the BAR result for the 
rest of the intervals as explained before. (Also note that Eq. ([1]) will yield AG instead of 
AF since we are dealing with NPT ensemble as mentioned at the end of sec. II). 

Figure [6] and [7] show the results for the BAR and the EBAR results using two simulations 
performed at A3 = 0.94 (state 0) and 1 (state 1). As seen in Fig. [6l the BAR result is in the 
small sample regime since J2o — Z^i 1- Also, we find that dAG/dC = — 1 at C = Cb 
from the slope of the plot in Fig. [TJ which again indicates a small sample regime.— As a 
result, the BAR result shows a large deviation from the actual value as seen in Fig. [71 On 
the other hand, the EBAR result (which corresponds to the maxima of \da'^/dC\ such that 
J2i > 1) is quite accurate as seen in Figs. [6l and [71 Note that as C approaches a large 
negative value C < —1170, the AG value approaches —33.8 ksT which corresponds to the 
FEP formula in state 0. Thus, the EBAR method yields more accurate results compared to 
both the single stage FEP method and the BAR method. 

In Fig. [HI we have tested the convergence of the results in the interval from A3 = 0.9 to 1. 
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(Note that the figure shows the total Gibbs free energy difference AGt for the entire path 
and all the results reported in the figure include a contribution of 60.75 ±3 ksT from Ai = 
to A3 = 0.9 computed by BAR method). As can be seen in the figure, BAR method requires 
12 simulations to obtain acceptable accuracy. On the other hand, EBAR (combined) method 
yields sufficiently accurate value of AGt with 4 simulations. Using the converged result in 
Fig. [8], we find that the contribution of the CM hysteresis to the total error in AGt is 
about ±2 ksT. We also found the thermodynamic integration using Gaussian Quadrature 
(GQ) method yields sufficiently accurate result with just 2 simulations in comparsion to 4 
simulations required by the EBAR method. In this case, GQ technique is effective because 
the integrand changes smoothly (although rapidly) as A3 approaches 1 and moreover GQ 
does not require evaluation of the integrand at A3 = 1 which is most prone to hysteresis. 

As for the convergence of the EBAR method, we note that for Ng = 2 (see Table [T] and 
Fig. [8]), the EBAR result (applied in the interval from A3 = 0.9 to 1) deviates significantly 
from the accurate value. This is because the external potential starts affecting the structure 
of the crystal phase significantly [see Eq. ([H])] for A3 < 0.9 (state 0) and hence the value of 
J2o becomes less accurate, as discussed in Sec. II. As we increase A^^, the EBAR (combined) 
result converges rapidly as seen in Fig. [HI This is because even for A3 = 0.99 the CM position 
is sufficiently localized and hence J2o evaluation is accurate. This can also be seen in the 
inset of Fig. [5l which shows that the CM hysteresis becomes appreciable only for A3 > 0.995. 

Finally, based on the value of AGt, we also computed the melting temperature (T^) by 
integrating the following equation at P* = 0:- 

' diAGT/ksT) ' 
dT 

J P,N 

where {■ ■ ■)l and (■ ■ ■)s denote the NPT ensemble averages for the liquid and crystal phases 
at the specified temperature. We found the value of T^ to be 1675 ± 5 K, based on EBAR 
(combined) result with Ng = 4. This is in close agreement with the value of 1678 K obtained 



kRT 



({U)T + P{V)L)-ms + P{V)s), (10) 



in Ref. 



I4J by crystal-melt coexistence simulations for Si(lOO) interface and the value of 



1691 ± 20 K obtained in Ref. 
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IV. SUMMARY 



In this work, we have shown that EBAR method efficiently calculates the free energy of 
the crystal phase due to an external potential without requiring use of a constraint on the 
translational degrees of freedom. In this method, we nullify the error incurred due to poor 
sampling of the perturbation energy in state (crystal phase without the external potential) 
by adjusting the value of the shift constant C [see Eq. ([1])] so that the error in estimated free 
energy difference is completely due to state 1 (the state in which an external potential acts on 
the crystal phase), where we expect the sampling of the perturbation energy to be sufficiently 
accurate. We have applied this technique to cleaving wall method, in which the crystal phase 
is subjected to a repulsive cleaving potential^ and confined fluid A-integration method,- in 
which the crystal phase is acted upon by attractive Gaussian potential wells located at the 
ideal crystal lattice sites. In both cases, we found that EBAR method yields accurate values 
with reasonable computational effort and and offers considerably improvement over both 
the single stage PEP method (in state 1) and the BAR method. Note that unlike the PEP 
method, the EBAR method utilizes information from both the ensembles. 

It must be stressed that the EBAR method is applicable only when the BAR result is in 
the small sample regime. Thus, the domains of applicability of the two methods are mutually 
exclusive. With regard to the convergence, we found that the EBAR method ceases to be 
accurate in the following limits: (i) When the external potential is too weak so that the CM 
position is not localized as seen in cleaving wall method (see Pig. H] and Table [T]) and (ii) 
when the external potential is too strong so as to affect the crystal structure significantly 
as in the case of constrained fluid A-integration method (see Pig. [Hand Table [T]). Both of 
these conditions increase the error the computation of J2i (corresponding to the state in 
which external potential is acting). The EBAR method is simple to implement since it only 
requires that the perturbation energies in the two states be sampled and does not depend 
upon the existence of a reversible path connecting the two states. 

In the constrained fluid A-integration method, thermodynamic integration by GQ tech- 
nique is found to be effective since the CM hysteresis is confined to the end of the integration 
path (see Pig. [5]). However, the GQ technique suffers from the inherent drawback of the 
TDI method in that it depends upon the existence of a reversible path between the two 
thermodynamic states. Also, testing the convergence of this technique is relatively expen- 
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sive since the abscissa values for higher number of integration points do not coincide with 
those corresponding to lower number of integration points.— These problems prevent general 
applicability of the GQ technique, as in the case of the cleaving wall method where the CM 
hysteresis occurs throughout the path (see Fig. [T]). 

We expect that EBAR method will also be useful in computing bulk crystal phase free 
energy by Einstein crystal method^ and in computing surface free energy of crystal phases^ 
where the CM hysteresis occurs. It seems possible to generalize the EBAR method to other 
free energy computations. Our initial calculations indicate that the EBAR method can also 
be applied to the cleaving of the crystal- melt interface (state 4 of the cleaving wall method), 
which is considered as a major challenge in the computation of the crystal-melt interfacial 
energy.-i^ Here the crystal-melt interface fluctuates in the absence of the external potential 
while it is held fixed when external cleaving potential is present.- It will be interesting to 
explore the applicability of EBAR method to the computation of the chemical potential 
by particle insertion-deletion technique.- Here the perturbation energies due to the particle 
insertion steps are sufficiently accurate while those due to the particle deletion steps are not 
sampled efficiently. This situation is similar to the problem considered in this article. 
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TABLE I: The number of simulations (Ng) and the free energy calculations obtained by the BAR 
(B), EBAR (E) and the combined (C) results. In the combined results {Ng > 2), the EBAR method 
is applied for the last interval while the BAR method is applied for the rest of the intervals. The 
second and the third column relate to the cleaving wall method, while the last two columns relate 
to the constrained fluid A-integration method. The data contained in this table is also plotted in 
Figs. Hand El 





z values 




AF/A (J/m2) 


A values 


AGr/kBT 


2 


0.75, 1.80 




0.041 ± 0.002 (E) 


0.9, 1.0 


28.6 ±7.1 (E) 








0.36 ±0.15 (B) 






4 


0.75, 0.80, 


0.85 


0.0455 ± 0.005 (C) 


0.9, 0.92, 0.94 


12.4 ±4.6 (C) 




1.80 






1.0 


-576 ± 193 (B) 


6 


0.75, 0.80, 


0.85 


0.0054 ± 0.009 (C) 


0.9, 0.92, 0.94 


9.0 ±4.4 (C) 




0.90, 1.00, 


1.80 


0.088 ± 0.024 (B) 


0.95, 0.96, 1.0 


-379 ± 130 (B) 


9 


0.75, 0.80, 


0.85 


0.0439 ± 0.004 (B) 


0.9, 0.92, 0.94 


5.5 ±4.4 (C) 




0.90, 1.00, 


1.10 




0.95, 0.96, 0.97 


-87.9 ± 34.4 (B) 




1.20, 1.50, 


1.80 




0.98, 0.99, 1.0 
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FIG. 1: The variation of the integrand in Eq. (5) as a function of z at T* = 0.0667, and pc = 
0.452(T~^ for cleaving of Si(lll) crystal phase. The inset shows the same plot with a fixed layer 
constraint. In applying this constraint the particles in the middle layers (layers numbering 42 and 
43) were immobilized. Note that due to the effect of the constraint on the free energy, the value of 
integrand at a given z is significantly larger in the inset plot. The error bars are seen when these 
are larger than the size of the symbols. 

FIG. 2: The crieteria for choosing values of C for the BAR and the EBAR methods. The data is 
generated from simulations done at z = 1.8 (state with the external potential) and z = 0.75 (state 
1 without the external potential). The BAR result corresponds to the condition that = X^i 
while the EBAR result corresponds to the maximum of \da'^/dC\ such that J2o is of order unity 
or higher. A local maximum also exists when J2i — 1; but is not seen in the figure because of its 
small magnitude. 

FIG. 3: The value of Free energy difference per unit area (in units of J/m?) obtained by Eq. ([1]) 
for the same set of data as that in Fig. [2j The inset shows more detailed comparison of the EBAR 
result with that of Ref. |8| (dashed line). From the slope of the plot, we get dAF/dC = —1 at 
C = Cb indicating a small sample regime for the BAR result.— At the two end points of the graph 
(C = — 10 and C = 220), the value of AF approach those obtained from the FEP formulas in the 
two ensembles. 

FIG. 4: Free energy change per unit area (in J/m^) resulting for cleaving of Si(lll) crystal phase 
without any constraint. The error bars are seen when these are larger than the size of the symbols. 
The abscissa represents the number of simulations {Ng) performed at different values of z (see 
Table U]). The horizontal line represents the result obtained using the BAR method in Ref. ^. The 
combined result (applicable for Ng > 2), represents a combination of the EBAR and the BAR 
methods as explained in the text. 

FIG. 5: The variation of the integrand in Eq. ^ as a function of A3 at P* = 0.0, T* = 0.0667, 
and N = 1000 for the forward and the backward runs for stage 3. The inset shows the region of 
maximum hysteresis. 
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FIG. 6: The crieteria for choosing appropriate values of C for the BAR and the EBAR methods. 
The data is generated from simulations performed at A3 = 0.94 (state 0) and A3 = 1.0 (state 1). 
The external Gaussian potential acts in state while it is absent in state 1 according to Eq. ([8]). 
The values of the Gibbs free energy difference AG between the two ensembles are plotted as a 
function of C in Fig. [71 The BAR result corresponds to the condition that J2o ~ Si while the 
EBAR result corresponds to the maximum of \da'^ /dC\ such that X^i — 1- 



FIG. 7: The value of Gibbs free energy difference (in units of ksT) between A3 = 0.94 and 1.0 
states. The result is obtained by Eq. ([1]) for the same set of data as that in Fig. [6l Note that the 
slope of the curve is — 1 at C = indicating a small sample regime for the BAR result. The 
inset shows more detailed comparison of the EBAR result with the accurate result (dashed line) 
obtained by inserting sevral intermediate steps between A3 = 0.94 and 1. At the two ends of the 
plot (C = —1170 and C = 0) AG values approach those corresponding to the FEP formulas in the 
two ensembles. 



FIG. 8: Total Gibbs Free energy difference (in units of fc^T) between the melt and the crystal 
phases obtained without applying any constraint. All of the results include a contribution of 
60.75 lb 3 k^T from Ai = to A3 = 0.9 as mentioned in the text. The abscissa represents the 
number of simulations performed at various values of A3 between 0.9 and 1. The horizontal line 
represents the converged value obtained by the BAR method with Ng = 12 (see the inset). The 
size of the error bar for this converged result is of the same magnitude as the combined result with 
Ns = 9. 
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